Method for automated window-level settings for magnetic resonance images

ABSTRACT

A method of automatically determining window-level settings in a digital image of a subject. The method includes the steps of: accessing the digital image; obtaining body part or coil type information; determining a region of interest (ROI) within the digital image; producing and processing a ROI-histogram; defining a W-L value for a linear or non-linear image intensity range re-mapping; and a linear or non-linear re-mapping the image intensity range.

FIELD OF THE INVENTION

[0001] The invention relates generally to the field of medical imaging,and in particular to digital image processing and display in MagneticResonance Imaging. More specifically, the invention relates to improvedimage representation for diagnostic purposes.

BACKGROUND OF THE INVENTION

[0002] Magnetic Resonance (MR) images are acquired during medicalexaminations on Magnetic Resonance Imaging (MRI) apparatus and stored indigital form in image archives. The digital data of the MR image can bestored and transferred to other image display, analysis, and archivingsystems by means of a computer network. The digital data is typicallytransmitted over a communication network using a DICOM Standard (DigitalImaging and Communications in Medicine (DICOM), Part 3: InformationObject Definitions, PS 3.3-2001).

[0003] The acquired pixel values in an MR image do not relate directlyto the brightness and contrast of display for diagnostic use of theimage, and frequently no information is provided about ideal windowwidth (contrast) and window level (brightness) that should be applied tothe image for this purpose. If the window width and level are present,very often this information is not optimal for diagnostic use. Thereforeradiologists or technicians usually need to manually adjust window widthand window level of these MR images (often referenced as window-level(W-L)) before the diagnosis can be performed on a display workstation.The manual adjustment procedures require the user's time and energy,particularly if a large number of MR slices are acquired for thepatient. Even a single study may include hundreds or thousands ofindividual images, grouped into several series for different spatialviews and/or acquisition parameters. Since the characteristics of theimages in a given series are similar, some medical personnel in order tobe efficient usually adjust W-L not for each slice individually but onlyfor the most clinically important one in each series (often this is aslice in the middle of the acquired volume, or close to it). However,when W-L values are set to the same value for all images in the sameseries, the results may not be optimal for all the images and mayrequire additional adjustment for diagnostic use.

[0004] Another consideration is that the range of gray scale used inconventional workstation monitors (typically, 8 bits) is oriented to thesensitivity of the human eye to small differences in light levels.However, the original image often has 1024 to 4096 or more gray levels(i.e., 10-12 or more bits). Therefore, a manual adjustment must be madewithout loosing diagnostically valuable features of the MR image. Thisis typically not a trivial task, particularly if there is more than onearea of diagnostic interest within the image, and also when those areassignificantly differ in intensity dynamic range. For these reasons,radiologists and MR technicians often differ in their subjectiveopinions about the W-L setting for the presentation of an image withbest image quality. For MRI, there is no unanimous agreement on the bestpresentation of this data, so a universal solution of the problem ofpicking W-L values does not exist. This is often referred to as the“Window-Level problem” or “W-L problem”.

[0005] Another consideration of the W-L problem are the characteristicsof MR images. For example: (1) MR images are measurements of very lowenergy quantum processes. As such, the images are often noisy, withmaximal and minimal pixel intensity values that vary from image toimage, and also within a particular image, e.g. spikes and backgroundnoise. (2) There are different types of possible MRI acquisitions, forexample, T1 or T2 weighted series. Each type of acquisition may requiredifferent parameters for optimal image display. (3) Image acquisitionmay be done with either (i) transmitting/receiving coils that surroundthe part of the patient being imaged (for example, head, extremity, orwhole body) or (ii) surface coils for structures such as the neck andspine. (4) The size and mass distribution of individual patients canaffect the values obtained from an MRI acquisition in complex ways. (5)Different users may prefer different “looks” of images with regard toimage brightness and contrast. (6) The display monitor type (e.g.,regular, flat panel), age (e.g., old, new), settings (e.g., gamma curve,etc.) have an effect on how an image area is displayed on a monitor.

[0006] One method to automatically set W-L values is to determine themaximal and the minimal intensity values of each MR image, and then tomap the values to an output range [0, 255] in a linear fashion. Thismethod, which is readily calculatable, does not work well for many MRimages, in particular, for those acquired with surface coils.

[0007] There have been prior art efforts to automate the manualadjustment procedures to address the Window-Level problem for differentmedical modalities, (e.g. X-ray, CT, MRI). Those efforts may beclassified into three groups. The first group comprises methods whereinthe operator participates in image quality evaluation to find the W-Lvalues. The second group comprises semi-automatic methods wherein theoperator provides some input or are based on values that have beendetermined empirically from operator experience. The third groupincludes adaptive and automatic methods.

[0008] With regard to the first group, one known approach is sometimesreferred to as “static contrast windowing”. This approach places themidpoint of the 8-bit dynamic range at the average intensity value ofthe region of interest (ROI) in the image, which is usually defined withoperator assistance. The intensity values that are below the contrastwindow are presented as black on the display and values above thecontrast window are displayed as white. Intensity values within thecontrast window are mapped to corresponding display brightness valuesusing a straight-line transfer curve. Contrast windowing is “static”when the window and the transfer curve are fixed for the conversion ofthe entire image data array. Details within the intensity range of thecontrast window are shown on the display using the 8 bits of dynamicrange. The remaining portion of the image, however, are displayed aseither “too white” or “too black”, with no details perceivable in thepresentation. When other areas of interest are not represented properlyat the selected W-L setting, more than one presentation of the image maybe needed, each with a different contrast window to properly representimportant anatomical details at all brightness levels.

[0009] Another approach sometimes used in X-ray CT practice employs atransfer curve having two contrast windows applied to a single image.This dual window approach uses two grayscale ramps or “windows” in asingle display, where one window is set to encompass the brightnesslevels around bone and the other window is set to cover the brightnesslevels around soft tissue. With this approach, the anatomy shown in CTimages is generally known to the user. As such, a trained user knows apriori what region of an image is soft tissue and what is bone, so thatthe simultaneous display of two gray ramps in the same image is notconfusing to the user. However, in any anatomic regions which areintermediate in brightness between the windowed soft tissue and bone,the X-ray CT image values may be displayed using parts of both grayscaleramps, which can be very confusing. This limits the applicability ofdual ramp presentation of CT. A common method of displaying such imagesuses separate presentations of the same image using the full bone rangeapplied to the image in one display area, and the soft tissue range isapplied to another image in a separate display area. However, neither adual-ramp presentation or multiple presentations with fixed W-L settingsmay be acceptable for MRI, since the digital values of images are notcalibrated to be independent of the individual examination as they arein CT. As discussed above, the proper Window-Level for any particularimage depends on a variety of factors that vary from image to image,from series to series, and from patient to patient.

[0010] One technique intended to enable the operator to select a staticwindow is described in U.S. Pat. No. 5,042,077 (Burke). The medicalimaging equipment produces an image of the subject under study andpresents a graph of the image's histogram, which indicates thedistribution of brightness levels of the image pixels. Using atrackball, the operator manipulates a contrast window which is displayedon the histogram and which enables the operator to select brightnessranges in the image for contrast enhancement.

[0011] With regard to the second group (i.e., the semi-automaticalgorithms), U.S. Pat. No. 5,900,732 (Felmlee) relates to an automaticwindowing method for MR images. The method includes producing ahistogram of the reconstructed image; removing histogram bins with lessthan a threshold (T) number of pixels; calculating a window level bydetermining the mean intensity of the remaining histogram bins; andcalculating the window width WW according to a fixed relationship. Whilethis method may have achieved a certain degree of success in itsparticular application, the method requires a manual step of setting thethreshold T, and therefore, requires the involvement of an MRItechnician or radiologist in the task of setting the W-L.

[0012] U.S. Pat. No. 5,835,618 (Fang) relates to a method for dynamicrange re-mapping for optimum image display. This method employs auser-defined parameter for both additive and multiplicative algorithms,as well as a “mixing” algorithm for the dynamic range re-mapping, thusrequiring an involvement of medical personnel into the technical detailsof adjusting W-L.

[0013] U.S. Pat. No. 5,357,549 (Maack) relates to a method of dynamicrange compression of an x-ray image wherein a non-linear mappingfunction is employed to determine the equalization value for each pixelin the image. This method's final output is dependent on a lowpassfiltered signal from the local neighborhood around each pixel. Adisadvantage of this method is the limitation of being able to displayonly one range of image data that is of primary interest with the bestgray scale presentation; a disadvantage particularly if more than onerange needs to be displayed simultaneously.

[0014] One non-linear method for image intensity correction in MRI isdescribed in Handbook of Medical Imaging, 2000. This method is based onintensity non-uniformity (INU) field estimation using thelog-transformed image histogram as a probability density function (PDF)approximation. This field is iteratively smoothed with a B-splinefunction and is used for intensity correction by means of mapping theoriginal histogram to the one sharpened by deconvolving a Gaussianblurring kernel.

[0015] With regard to the third group (i.e., adaptive or automaticalgorithms), one known method for mapping data brightness levels to adisplay having a limited dynamic range is known in the art as adaptivecontrast enhancement. One variant is “adaptive histogram equalization”(AHE). Unlike methods which are “static”, the AHE method does not employa fixed contrast window for the entire image. Instead, the AHE methodlooks at each datum intensity value in the acquired data array one at atime and compares it with the values in a local surrounding spatialarea, or “context region”. The length and width of the context regionmay, for example, range from one sixth to one sixtieth of the length andwidth of the entire image data array. While there are many variations onthe calculations employed with this method, generally the centered datumvalue is mapped to a display brightness, which provides contrast withrespect to the other data values that are generally within the samecontext region. The calculations are performed, in principle, at eachpixel location in the image data array with respect to its surroundingcontext region and the method is, therefore, computationally intensive.The AHE method and some variations are described in “Algorithms ForAdaptive Histogram Equalization”, Pizer, S. M., Austin, J. D. et al.,SPIE Vol. 671, Physics and Engineering of Computer MultidimensionalImaging and Processing, 1986.

[0016] U.S. Pat. No. 5,305,204 (Ohhashi) relates to a digital imagedisplay apparatus with automatic window level and window widthadjustment. With this apparatus, a pixel value of a digital image datais converted into brightness by using a suitable conversion function,and displaying an image within an optimum window. This pure neuralnetwork based method may not be reliable or accurate when encountering anew kind of input image which has not been seen in the training set ofthe neural network.

[0017] Another method which employs a neural networks for automaticwindow-level adjustment in MRI is disclosed in U.S. Pat. No. 5,995,644(Lai). A further method which employs a neural networks is disclosed inU.S. Pat. No. 6,175,643 (Lai).

[0018] U.S. Pat. No. 5,268,967 (Jang) relates to a method for automaticforeground and background detection in digital radiographic images. Themethod includes the step of detecting edges in the image signal based ona morphological edge detector.

[0019] Accordingly, while the apparatus and methods described above mayhave achieved certain degrees of success in their particularapplications, the manual methods in W-L adjustment can lead tovariations in image quality because individual operators have differentopinions regarding a preferred presentation of specific MRI images data.In addition, the manual adjustment requires time, therefore reducingthroughput of the medical imaging system for diagnosis. Thesemi-automated methods employ manual selection of parameters or requirethe operator to evaluate one or more criteria of an image or study.Prior art adaptive methods may also need some operator participation orat least some prior information to train the system, which producesfinal W-L parameters for display.

[0020] As such, there exists a need for an apparatus and method whichovercomes the problems described above. The present invention isdirected to overcoming the problems described above. In particular, thepresent invention is directed to a fully automated method. The method isbased on modeling radiologist's criteria for the appearance of MR imagesand comprises of a sequence of steps in W-L adjustment. This methodperforms both spatial image and histogram analysis and processing, andcan include the steps of obtaining of the body part or coil typeinformation, extraction of the region of interest (ROI) within theoriginal digital image, a W-L definition based on pixel valuedistribution analysis within extracted ROI, and linear or non-linearre-mapping of the original image range of intensities to the gray scalerange of the display.

SUMMARY OF THE INVENTION

[0021] An object of the present invention is to provide a fullyautomated method for window-level setting for MR images.

[0022] Another object of the present invention is to provide such afully automated method which requires no or minimal adaptation ortraining for the method itself or for the users.

[0023] A further object of the present invention is to provide such afully automated method, which provides good image quality.

[0024] These objects are given only by way of illustrative example, andsuch objects may be exemplary of one or more embodiments of theinvention. Other desirable objectives and advantages inherently achievedby the disclosed invention may occur or become apparent to those skilledin the art. The invention is defined by the appended claims.

[0025] The method of the present invention analyzes a comprehensive setof image characteristics collected in spatial, spectral and histogramdomains and header information. In addition, the method of the presentinvention includes two selectable options (linear or non-linear) topresent images depending on the radiologist's personal perception anddiagnosis needs. The non-linear window-level method of the presentinvention provides different image areas within ROI that might beclinically valuable to be displayed at a time, which might be impossibleusing only the linear method. The method is fully automated and does notneed any adaptation or training for the method itself or for the users.

[0026] According to one aspect of the invention, there is provided amethod of automatically determining window-level settings in a digitalimage of a subject. The method comprises the steps of: accessing thedigital image; obtaining body part or coil type information; determininga region of interest (ROI) within the digital image; producing andprocessing a ROI-histogram; defining a W-L value for a linear ornon-linear image intensity range re-mapping; and a linear or non-linearre-mapping the image intensity range.

BRIEF DESCRIPTION OF THE DRAWINGS

[0027] The foregoing and other objects, features, and advantages of theinvention will be apparent from the following more particulardescription of the preferred embodiments of the invention, asillustrated in the accompanying drawings.

[0028]FIG. 1 shows a block diagram illustrating a process flow inaccordance with the present invention for MR image input, processing,and display or archive, including an automatic window-level method.

[0029]FIG. 2 shows a block diagram illustrating an algorithm inaccordance with the present invention for determining a coil type usedin acquiring an MR image.

[0030]FIGS. 3a-3 d illustrate an image row and column mean value vectorsrelationship in FIG. 3a for head sagittal view (surrounding coil); inFIG. 3b for head coronal view (surrounding coil), in FIG. 3c for spinesagittal view (surface coil); and in FIG. 3d for spine axial view(surface coil).

[0031]FIG. 4 shows a block diagram illustrating an image ROI extractionalgorithm in accordance with the present invention.

[0032]FIG. 5 shows a block diagram illustrating an algorithm forsmoothed ROI-histogram analysis and window-level values calculation forlinear and non-linear intensity re-mapping in accordance with thepresent invention.

[0033]FIG. 6 shows a block diagram illustrating a non-linear option forintensity re-mapping in accordance with the present invention.

[0034]FIGS. 7a-7 d illustrate histograms and re-mapping functionscalculated for a sagittal head image in FIG. 7a the original imagehistogram; in FIG. 7b the smoothed ROI-histogram (original histogram isshown in dotted line); in FIG. 7c the smoothed ROI-histogram thresholdedfor non-linear intensity re-mapping method (threshold for linear methodis shown in dotted line); and in FIG. 7d the re-mapping transferfunctions.

[0035]FIGS. 8a-8 d illustrate histograms and re-mapping functionscalculated for a coronal head image within in FIG. 8a the original imagehistogram; in FIG. 8b the smoothed ROI-histogram (original histogram isshown in dotted line); in FIG. 8c the smoothed ROI-histogram thresholdedfor non-linear intensity re-mapping method (threshold for linear methodis shown in dotted line); and in FIG. 8d the re-mapping transferfunctions.

[0036]FIGS. 9a-9 d illustrate histograms and re-mapping functionscalculated for an axial head image within in FIG. 9a the original imagehistogram; in FIG. 9b the smoothed ROI-histogram (original histogram isshown in dotted line); in FIG. 9c the smoothed ROI-histogram thresholdedfor non-linear intensity re-mapping method (threshold for linear methodis shown in dotted line); and in FIG. 9d the re-mapping transferfunctions.

[0037]FIGS. 10a-10 d illustrate histograms and re-mapping functionscalculated for sagittal spine image within in FIG. 10a the originalimage histogram; in FIG. 10b the smoothed ROI-histogram (originalhistogram is shown in dotted line); in FIG. 10c the smoothedROI-histogram thresholded for non-linear intensity re-mapping method(threshold for linear method is shown in dotted line); and in FIG. 10dthe re-mapping transfer functions.

[0038]FIGS. 11a-11 d illustrate histograms and re-mapping functionscalculated for axial spine image within in FIG. 11a the original imagehistogram; in FIG. 11b the smoothed ROI-histogram (original histogram isshown in dotted line); in FIG. 11c the smoothed ROI-histogramthresholded for non-linear intensity re-mapping method (threshold forlinear method is shown in dotted line); and in FIG. 11d the re-mappingtransfer functions.

[0039]FIGS. 12a-12 b illustrate an output MR sagittal head image withnon-linear re-mapped intensities (FIG. 12a) and corresponding ROI binarymask (FIG. 12b).

[0040]FIGS. 13a-13 b illustrate an output MR coronal head image withnon-linear re-mapped intensities (FIG. 13a) and corresponding ROI binarymask (FIG. 13b).

[0041]FIGS. 14a-14 b illustrate an output MR axial head image withnon-linear re-mapped intensities (FIG. 14a) and corresponding ROI binarymask (FIG. 14b).

[0042]FIGS. 15a-15 b illustrate an output MR sagittal spine image withnon-linear re-mapped intensities (FIG. 15a) and corresponding ROI binarymask (FIG. 15b).

[0043]FIGS. 16a-16 b illustrate an output MR axial spine image withnon-linear re-mapped intensities (FIG. 16a) and corresponding ROI binarymask (FIG. 16b).

DETAILED DESCRIPTION OF THE INVENTION

[0044] The following is a detailed description of the preferredembodiments of the invention, reference being made to the drawings inwhich the same reference numerals identify the same elements ofstructure in each of the several figures.

[0045] The present invention recognizes that the spatial intensitydistribution of an image plays an important role in the adjustment ofthe display window parameters, and different organs may have similarhistograms but might be windowed differently. The present invention alsorecognizes the need to focus on the analysis of radiologist/operatoractions while adjusting W-L.

[0046] The present invention provides a method for automatic calculationof the window width and window level values and intensity rangere-mapping for MR images for display and/or archive. The algorithmincludes five steps: body part or coil type information extraction, ROIextraction, processing of histogram calculated within the ROI,window-level values definition, and linear or non-linear originalintensity re-mapping for image display. Body part information might beobtained from image header or might be estimated by coil type detectionfrom the image itself, if header information is not available. ROIextraction comprises thresholding of image pixel values, then groupingthe discriminated values by means of their external bounding followed bymorphological filtering of the ROI boundary. ROI-histogram is calculatedwithin the extracted ROI and the histogram is then smoothed. Windowwidth and level are defined by means of body part/coil type dependentthresholding and adjusting of the processed ROI-histogram for furtherboth linear and non-linear image intensity re-mapping methods. Linearre-mapping is well known standard method and is done by simpleprojection of image pixel values within the defined window to a displaygray scale range. Non-linear re-mapping presumes that image pixels aretransformed using a non-linear weighting function. In both re-mappingmethods, the image pixel values that are less/greater than window widthlower/upper limit are set to black/white respectively. The non-linearmethod can be used as an alternative user selectable option and can beuseful when multiple areas with different dynamic ranges within the ROIare needed to be displayed, especially for images with a very broadoriginal range (e.g., more than 12 bits).

[0047] Applicants have recognized two goals of a radiologist/technicianin adjusting image brightness and contrast: first, within clinicallyimportant areas to reach the right balance between the amount of darkand bright spots (using those limits as landmarks), and second, todisplay the rest of the chosen area with maximum contrast withinavailable gray scale. The clinically important area is often a globalROI (i.e., the whole body part), but sometimes it might be a sub-regionof the image (i.e., a specific area within a body part). In the lattercase, particularly when different local areas have different intensityranges, it might be necessary to produce several images with differentW-L values or to find a compromise setting to display them together,depending on body part.

[0048] The present invention employs several steps of image analysis andprocessing: body part or coil type information obtaining, global ROIextraction, intensity distribution analysis within ROI, W-L valuescalculation, image intensity range re-mapping and image display/archive.

[0049] Referring to FIG. 1, there is shown a block diagram illustratinga process flow in accordance with the present invention for MR imageinput, processing, and display or archive, including an automaticwindow-level method. More particularly, a MR digital image is acquired(block 100) and the MR image is transmitted to a workstation memory(block 102). At block 104, a body part or coil type information isobtained. Then, the image ROI is extracted (block 106). A histogramwithin the ROI is then calculated and smoothed (block 108). Window-levelvalues are defined for both the linear and non-linear methods (block110). If a linear method is determined (block 112), then a linear imageintensity range re-mapping occurs (block 114). If a non-linear method isdetermined, then a non-linear image intensity range re-mapping occurs(block 116). The linear method might be used as a default while thenon-linear method may be preferable in case of multiple local ROI orwhen a very broad original range of intensities exists in the image.Finally, the image can be displayed or/and archived, or the calculatedW-L values may be stored in the image header for transmission to otherstorage or display systems (block 118).

[0050] With regard to blocks 100 and 102, the acquisition of an MRdigital image and its transmission to a workstation memory are wellknown to those skilled in the art.

[0051] Block 104 comprises the extraction of body part or coil typeinformation. It is noted that this information may or may not beexplicit in the additional data that describes the image. That is,although DICOM includes tags that can describe the body region andacquisition conditions, some of these tags are optional, may be empty,or may be filled inconsistently. For example, DICOM Tags Receiving Coil(0018,1250) & Transmitting Coil (0018,1251) are both DICOM Type 3, andas such, they are completely optional, or if present may be empty. Inaddition, there are no defined terms for these values, so it is notpossible to use them to determine consistently good information on theacquisition details. But if the information in those tags is availablethen it can be used for input parameters (i.e., body part and view).

[0052] Applicants have noted that, with many MRI studies, images of thecervical, thoracic, and lumbar spine regions require a different W-Lcalculation than images from other parts of the body, including thehead, chest, abdomen, and extremities. Applicants believe that adifference is that spine images are acquired using a surface coilwhereas other body regions are taken with coils that surround the bodyregion, such as a head coil.

[0053] Images taken using surrounding coils provide consistentbrightness for a given tissue throughout the image area. For example,the average brightness of similar brain tissue is consistent over animage acquired with a head coil. An analysis of the histogram of such animage is simplified by the uniformity of values in the image for thesame tissue through the image region.

[0054] An effect of a surface coil is that the anatomy close to the bodysurface is relatively bright, and parts of the image farther from thecoil are dark. The brightness of a tissue depends not only on thecharacteristics of that tissue but also on the distance relative to thecoil. This results in bright regions in the image generally found in aregion under the skin and parallel to the anatomy of interest. Thereforethe presentation of such an image to show the tissues of interestdepends on knowledge of how the coil is positioned relative to the bodyregion to be examined.

[0055] Other specialized coils may be employed for some kinds ofimaging, including breast imaging and also MRI of the jaw fortemporal-mandibular joint (TMJ) evaluations. These also producelocalized regions in the images that are especially bright near the coilpositions.

[0056] The present invention employs an algorithmic method thatdetermines if it is likely that the MR image was acquired using asurface coil or with a surrounding coil. This method is based on imagesymmetry analysis. These data are then used as input to the analysismethod that determines the actual W-L values for the image if there wasno information in the image header. In addition, it is generally thecase that all images in single series of an MRI examination are acquiredwith the same physical configuration. Thus, even if the determination ofsurfaces vs. surrounding coil is not accurate for each individual image,applying the algorithm to each image in the series will provide a strongindication of the type of coil used for the entire series.

[0057] If a symmetrical image is considered as a set of random variablesobtained by sampling a stationary stochastic process (MRI signalacquisition), then its statistical properties should not depend greatlyon the direction of sampling the image's values. To verify this, anestimation can be made of the variation of local mean values, calculatedfor each image profile individually, from global one (for the wholeimage) depending on direction. The algorithm for coil detection isdependent on this characteristic.

[0058] Referring now to FIGS. 2, for a particular MR image, the coildetection algorithm of block 104 is now discussed. More particularly,block 104 is comprised of the steps disclosed by blocks 200-210.

[0059] With the input image in the workstation memory (block 102), foreach row of the image, a mean value is computed to obtain a vector ofrow mean values and then find that row vector maximal deviation (RMD)from an image global mean value (GMV) (block 200). Then, at block 202,for each column of the image, a mean value is computed to obtain avector of column mean values and then find that column vector maximaldeviation (CMD) from the GMV.

[0060] If the CMD is greater than a predetermined value (block 204),then the data is determined to a surface coil image (block 206), and theimage ROI extraction is accomplished (block 106). This predeterminedvalue is preferably an empirically determined constant times the RMD,more particularly, 2 times RMD.

[0061] If the RMD is greater than a predetermined value (block 208),then the data is determined to a surface coil image (block 206), and theimage ROI extraction is accomplished (block 106). This predeterminedvalue is preferably an empirically determined constant times the CMD,more particularly, 2 times CMD.

[0062] If neither the CMD nor the RMD is greater than the predeterminedvalue (i.e., blocks 204 and 208 are not met), then the data isinterpreted as an image acquired with a surrounding coil (block 210).

[0063] It is noted that that the empirically determined constant ofabout 2 for image asymmetry factor in blocks 204 and 208 has been foundto give acceptable results.

[0064] Accordingly, the output of the step of block 104 is either (i)body part (and/or view) obtained from the MR image header or (ii) coiltype information extracted from image analysis. In a case of acorrelation of both outputs, the header information is preferable,however in case of a conflict, the coil type information has higherpriority. For example, if the image header information was incorrectlycompleted.

[0065]FIGS. 3a-3 d illustrate an image row and column mean value vectorsrelationship, i.e. those vectors deviation from GMV. FIG. 3a illustratesa head sagittal view (surrounding coil). FIG. 3b illustrates a headcoronal view (surrounding coil). FIG. 3c illustrates a spine sagittalview (surface coil). FIG. 3d illustrates a spine axial view (surfacecoil). The deviations of mean value vectors from GMV for head images donot differ much, illustrating image symmetry, while large difference inmean value vector deviations for spine images indicates image asymmetry.

[0066] The next step in the method of the present invention is globalROI extraction. Accordingly, referring now to FIG. 4, the image ROIextraction block 106 is now discussed. More particularly, block 106 iscomprised of the steps disclosed by blocks 302-310.

[0067] As previously considered, image profiles are random variables.Since an image profile is a 1-D spatial distribution of intensities(unsigned function), its frequency spectrum has a peak at a frequencyequal to DC. The 2-D frequency spectrum of such image will also have apeak at DC, which indicates a presence of constant shift of all imagevalues.

[0068] It is known from statistics (refer to Fuller W. A. “Introductionto Statistical Time Series”, Wiley, John & Sons, 1995.) and signalprocessing (refer to Wahl F. M. “Digital image signal processing”,Artech, 1987.) theory that subtracting the random variable mathematicalexpectation (mean value) from all its values will remove such shift.This technique is acceptable for signed functions (signals) to improvethe signal-to-noise ratio (SNR) or to calculate other statisticalcharacteristics. However, this technique cannot be used the same way forunsigned (image intensity) values/functions. Moreover, image pixelvalues less than the mentioned above mean value correspond to imagebackground and very low intensity values that might be of some interest,or, that is, they correspond to image histogram bins of very low or zerointensity values. If those bins have very high values, they can make theimage histogram less informative (as will be shown below with referenceto FIGS. 7a and 11 a) and sometimes even useless for further analysis orprocessing (as will be shown below with reference to FIGS. 8a, 9 a, and10 a).

[0069] Some prior art methods described above that remove very lowintensity values in the image histogram seem to be arbitrary and notbased on the expected characteristics of MRI data in spatial orfrequency domain. At block 104, the body part or coil type informationis obtained, and then an image GMV is determined (block 302). Moreparticularly, in the present invention, instead of mean valuesubtraction, the GMV is employed as a threshold in the spatial domain todifferentiate pixel values greater than it from all other pixels (block304) for further ROI extraction and its histogram analysis. In anotherwords, using such thresholding (segmentation) techniques for ROIextraction in space domain, the SNR is improved for the informative partof image histogram, thereby eliminating DC shift in the image spectrum.This image segmentation is performed in binary form, i.e. pixels abovethe threshold are set to black and background pixels are set to white toobtain an initial binary mask. The extracted black areas/pixelscorrespond to low-medium, medium and high intensities and are presumedto be the areas of clinical interest that would be expected from thecharacteristics of typical MRI modality signal-to-noise ratio. Thoseareas or pixels of interest (POI) are grouped by means of their externalbounding (like convex shell) (block 306). The area inside the externalboundary (EB) is set to black followed by standard (see Gonzales R. C.,Woods R. E. “Digital Image Processing”, Addison Wesley Longman, 2001.)morphological filtering (MF) to exclude spikes and the influence ofbackground noise and also to include some areas that are possiblyspatially separated from the main area of the segmented image (block308), i.e., the islands are joined to the mainland after flooding. Aglobal ROI is then obtained as a binary mask, (as will be shown belowwith reference to FIGS. 11b, 12 b, 13 b, 14 b, 15 b, and 16 b) formed bypixels of interest external boundary morphologically filtered (POI EBMF, block 310).

[0070] The next step in the method of the present invention isROI-histogram calculation and smoothing. Accordingly, block 108 is nowdiscussed (block 108 is shown in FIGS. 1, 4, and 5). The histogram(ROI-H) is calculated from the input image within the extracted ROI andsmoothed. The ROI-H histogram needs to be smoothed to avoid spikes thatwould affect further analysis. This smoothing might be done withwell-known methods such as standard (see Gonzales, Woods, “Digital ImageProcessing”, 2001) morphological or Gaussian or lowpass filtering. Ifthe lowpass filtered ROI-H-S histograms (as will be shown below withreference to FIGS. 7b, 8 b, 9 b, 10 b, and 11 b) are compared with theoriginal histograms (as will be shown below with reference to FIGS. 7a,8 a, 9 a, 10 a, and 11 a), it is observed that the large valuescorresponding to image background and low intensities are dramaticallychanged. These values did not disappear completely though, because theextracted ROI might also contain some low intensity values.

[0071] The next step in the method of the present invention is definingW-L values for further steps of linear and non-linear intensity rangere-mapping (block 110). Accordingly, referring now to FIG. 5, block 110includes the steps disclosed by blocks 402-410.

[0072] As previously described, when radiologists/technicians adjustW-L, they are balancing the amount of black and white values with thenumber of the remaining pixel data, trying to display as many imagedetails as possible within available gray scale. The ROI-H-S histogramsare more informative than the original ones but they may still includevery “low populated” (i.e., low informative) regions that may be removedto determine the most valuable intensity range.

[0073] Accordingly, to balance the amount of black and white pixels withthe rest pixel data the present invention employs smoothed ROI-histogramthresholding, which is body part/view dependent, to define W-L valuesfor further linear and non-linear re-mapping of original intensities.The intensity range differences for linear and non-linear re-mappingmethods exist, because the non-linear method can display a broaderintensity range than the linear method without a big decrease in theimage contrast.

[0074] As described above, surface-coil images (usually used for spineimaging) tend to have a larger dynamic range than the body (surrounding)coil images.

[0075] For the non-linear re-mapping method, the present inventionemploys a median threshold for all body parts in the case of noinformation about the view presented in the image (block 402). However,if such information is available, then for spine axial data view thethreshold is chosen to be the mean value of smoothed ROI-histogram, andfor coronal view head images there is no threshold. After thresholding,the histogram may need some additional adjustment in some regions. Theremight be several histogram portions that exceed the threshold: primary(the widest one), secondary (more narrow), etc., which may be separatedby gaps. If the intensity gap between the primary and any neighboringhistogram segment is narrower than the width of the neighboring, thenthis segment is attached to the primary one. Such a histogram adjustment(block 404) prevents the loss of any valuable intensity range fornon-linear intensity re-mapping.

[0076] For the linear re-mapping method, the present invention employs amean value threshold for spine (or surface-coil) images and median valuethreshold for images of other body parts (block 408) to extract the mostinformative portion of the ROI-H-S histogram (as will be shown belowwith reference to FIGS. 7c, 8 c, 9 c, 10 c, and 11 c, dotted line).

[0077] The intensity range between the first and the last values ofthresholded and adjusted ROI-H-S-NLT-A (block 406) and ROI-H-S-LT (block410) histograms is considered as the window width and the middle pointis considered as the window level. Those W-L values are defined forlinear and non-linear re-mapping of the image intensity range. It mightbe useful to store the linear method values since a linear approach isstandard practice and is widely used.

[0078] A last step for the present invention is a linear or non-linearre-mapping of the image intensity range. For both methods, the imagepixels values (1) that are less/greater than a lower/upper limit ofrespective window are windowed (I^(w)), i.e., they are set toblack/white respectively.

[0079] A linear method of image range re-mapping to display gray scalerange is a standard and known linear transform of image pixel valueswithin the defined window width (as will be shown with reference toFIGS. 7d, 8 d, 9 d, 10 d, 11 d, dotted line).

[0080] A non-linear re-mapping option can be used when there are severalareas of interest with different intensity ranges within the defined ROIor when the image needs to be displayed with more contrast for detailedanalysis (referring again to FIG. 1, block 116).

[0081] A proposed non-linear option for the present invention is animage intensity range re-mapping based on automated adjustment of lightflux going through the ROI area, which might be reformulated in terms ofa known histogram equalization technique (refer to Gonzales and Woods,Digital Image Processing, 2001). This technique is applied to ROI or,which is equivalent, to the histogram window defined in a previous stepfor the non-linear method. This non-linear method is implemented asillustrated in FIG. 6, and includes the steps shown in blocks 502-506.

[0082] At block 502, the cumulative distribution function (CDF) iscalculated and normalized for the processed ROI-H-S-NLT-A histogramwithin the defined window. This CDF is then used as a non-lineartransfer function for intensity re-mapping (as will be shown withreference to FIGS. 7d, 8 d, 9 d, 10 d, and 11 d, solid line). For thispurpose weighting function coefficients W_(i) are calculated (block 504)according to the relationship:

W=1+CDF _(i)−(I _(i) −I _(L))/Δ  (Equation 1)

[0083] where i is an integer index within the interval equal to thedefined window width, CD_(i) is a normalized CDF value, I_(i) is theimage pixel value, I_(L) is the lower limit of the defined window, andΔl is the defined window width.

[0084] Next, all windowed image pixels I^(w) are weighted using thefunction W(I) to get final image pixel values I′ (block 506) fordisplay:

I′=W(I)·I ^(w)  (Equation 2)

[0085] $I_{i}^{w} = \{ \begin{matrix}I_{L} & {{i \leq L} = \text{window lower limit}} \\{I_{i},} & {L < i < U} \\{I_{U},} & {{i \geq U} = \text{window upper limit}}\end{matrix} $

[0086] where I^(W) is windowed image, i.e.

[0087] Finally, the linearly or non-linearly intensity re-mapped imageis displayed at the workstation for diagnosis, stored, or archived withthe W-L values calculated for linear intensity re-mapping method to beincluded with the image in the archive (block 118). Usually, imagespresented with linear re-mapped intensities appears to be a bit brighterin general than non-linear W-Led one. This provides flexibility to theproposed method so that a user's personal preferences can be used toselect the desired presentation.

[0088]FIGS. 7-11 more particularly illustrate the method of the presentinvention.

[0089]FIGS. 7a-7 d illustrate histograms and re-mapping functionscalculated for a sagittal head image. FIG. 7a shows the original imagehistogram. FIG. 7b shows the smoothed ROI-histogram (the originalhistogram is shown in dotted line). FIG. 7c shows the smoothedROI-histogram thresholded for non-linear intensity re-mapping method(threshold for linear method is shown in dotted line). FIG. 7d shows there-mapping transfer functions.

[0090]FIGS. 8a-8 d illustrate histograms and re-mapping functionscalculated for a coronal head image. FIG. 8a shows the original imagehistogram. FIG. 8b shows the smoothed ROI-histogram (original histogramis shown in dotted line). FIG. 8c shows the smoothed ROI-histogramthresholded for non-linear intensity re-mapping method (threshold forlinear method is shown in dotted line). FIG. 8d shows the re-mappingtransfer functions.

[0091]FIGS. 9a-9 d illustrate histograms and re-mapping functionscalculated for an axial head image. FIG. 9a shows the original imagehistogram. FIG. 9b shows the smoothed ROI-histogram (original histogramis shown in dotted line). FIG. 9c shows the smoothed ROI-histogramthresholded for non-linear intensity re-mapping method (threshold forlinear method is shown in dotted line). FIG. 9d shows the re-mappingtransfer functions.

[0092]FIGS. 10a-10 d illustrate histograms and re-mapping functionscalculated for a sagittal spine image. FIG. 10a shows the original imagehistogram. FIG. 10b shows the smoothed ROI-histogram (original histogramis shown in dotted line). FIG. 10c shows the smoothed ROI-histogramthresholded for non-linear intensity re-mapping method (threshold forlinear method is shown in dotted line). FIG. 10d shows the re-mappingtransfer functions.

[0093]FIGS. 11a-11 d illustrate histograms and re-mapping functionscalculated for an axial spine image. FIG. 11a shows the original imagehistogram. FIG. 11b shows the smoothed ROI-histogram (original histogramis shown in dotted line). FIG. 11c shows the smoothed ROI-histogramthresholded for non-linear intensity re-mapping method (threshold forlinear method is shown in dotted line). FIG. 11d shows the re-mappingtransfer functions.

[0094]FIGS. 12a-12b illustrate an output MR sagittal head image withnon-linear re-mapped intensities (FIG. 12a) and corresponding ROI binarymask (FIGS. 12b).

[0095]FIGS. 13a-13 b illustrate an output MR coronal head image withnon-linear re-mapped intensities (FIG. 13a) and corresponding ROI binarymask (FIG. 13b).

[0096]FIGS. 14a-14 b illustrate an output MR axial head image withnon-linear re-mapped intensities (FIG. 14a) and corresponding ROI binarymask (FIG. 14b).

[0097]FIGS. 15a-15 b illustrate an output MR sagittal spine image withnon-linear re-mapped intensities (FIG. 15a) and corresponding ROI binarymask (FIG. 15b).

[0098]FIGS. 16a-16 b illustrate an output MR axial spine image withnon-linear re-mapped intensities (FIG. 16a) and corresponding ROI binarymask (FIG. 16b).

[0099] The method of the present invention has been implemented, and asimplemented, performance in the range of about 0.05 sec (on a PentiumIII/500 MHz PC) has been accomplished to process each MRI slice of size512×512 pixels.

[0100] A computer program product may include one or more storagemedium, for example; magnetic storage media such as magnetic disk (suchas a floppy disk) or magnetic tape; optical storage media such asoptical disk, optical tape, or machine readable bar code; solid-stateelectronic storage devices such as random access memory (RAM), orread-only memory (ROM); or any other physical device or media employedto store a computer program having instructions for controlling one ormore computers to practice the method according to the presentinvention.

[0101] The invention has been described in detail with particularreference to a presently preferred embodiment, but it will be understoodthat variations and modifications can be effected within the spirit andscope of the invention. The presently disclosed embodiments aretherefore considered in all respects to be illustrative and notrestrictive. The scope of the invention is indicated by the appendedclaims, and all changes that come within the meaning and range ofequivalents thereof are intended to be embraced therein.

What is claimed is:
 1. A method of automatically determiningwindow-level settings in a digital image of a subject, comprising thesteps of: accessing the digital image; obtaining body part or coil typeinformation; determining a region of interest (ROI) within the digitalimage; producing and processing a ROI-histogram; defining a W-L valuefor a linear or non-linear image intensity range re-mapping; and linearor non-linear re-mapping the image intensity range.
 2. The method ofclaim 1, wherein the step of obtaining coil type information comprisesimage symmetry analysis.
 3. The method of claim 2, wherein the imagesymmetry analysis comprises the steps of: calculating a vector of rowmean values; determining a maximal deviation of the row vector from animage global mean value; calculating a vector of column mean values;determining a maximal deviation of the column vector from the imageglobal mean value; and employing a predetermined value image asymmetryfactor.
 4. The method of claim 3, wherein the predetermined value imageasymmetry factor has a value of about
 2. 5. The method of claim 1,wherein the step of determining a region of interest (ROI) comprisesimage segmentation in a spatial domain.
 6. The method of claim 5,wherein the image segmentation comprises the steps of: determiningpixels of interest (POI) employing an image global mean value (GMV) as athreshold; grouping the POI by defining an external boundary (EB); andmorphologically filtering the region within the EB to define the ROI. 7.The method of claim 1, wherein the step of producing and processing aROI-histogram comprises the steps of: calculating a histogram within theROI; and smoothing the ROI-histogram with a lowpass, morphological orGaussian filter.
 8. The method of claim 1, wherein the step of defininga W-L value comprises the step of thresholding and adjusting theROI-histogram.
 9. The method of claim 8, wherein the step ofthresholding is accomplished by employing a threshold having a valuedifferent for the linear image intensity range re-mapping than for thenon-linear image intensity range re-mapping.
 10. The method of claim 9,wherein the step of adjusting is accomplished by the steps of attachingvaluable portions of the thresholded ROI-histogram for the non-linearimage intensity range re-mapping to its widest portion when the valuableportions are separated from the widest one by gaps in intensity.
 11. Themethod of claim 1, wherein the step of defining a W-L value for anon-linear image intensity range re-mapping comprises the steps of:smoothing, thresholding, and adjusting the ROI-histogram; defining awindow width as the difference between a last and a first non-zero valueof the smoothed, thresholded and adjusted ROI-histogram; and defining awindow level as a middle point between the last and first non-zerovalues of the smoothed, thresholded and adjusted ROI-histogram.
 12. Themethod of claim 1, wherein the step of non-linear re-mapping the imageintensity range comprises the step of windowing the input image pixelsby setting all input image pixel values to black/white for image pixelsthat are less/greater than a predetermined lower/upper window limit,respectively.
 13. The method of claim 1, wherein the step of non-linearre-mapping comprises the steps of: calculating a non-linear intensitytransfer function; calculating non-linear weighting functioncoefficients; and non-linear transforming of input image pixels.
 14. Themethod of claim 13, wherein the step of calculating a non-linearintensity transfer function is based on a cumulative distributionfunction (CDF) calculation from a smoothed, thresholded and adjustedROI-histogram within a defined window for non-linear image intensityrange re-mapping.
 15. The method of claim 13, wherein the non-linearweighting function coefficients W_(I) are calculated by therelationship: W _(i)=1+CDF _(i)−(I _(i−I) _(L))/ΔI where i is aninteger, CDF_(i) is a normalized CDF value, I_(i) is an image pixelvalue, I_(L) is a window lower limit, and ΔI is a defined window width.16. The method of claim 13, wherein the step of non-linear transformingis accomplished according to the relationship: I′=W(I)·I ^(w) whereI′are transformed image pixel values, W(I) is a weighting function,I^(w) are windowed input image pixel values.
 17. The method of claim 1,wherein the digital image is a magnetic resonance image acquired using amagnetic resonance imaging (MRI) system.
 18. The method of claim 1,further comprising the step of displaying, archiving, or storing thedigital image.
 19. The method of claim 1, further comprising the step ofdisplaying, archiving, or storing the W-L value.
 20. A computer storageproduct having at least one computer storage medium having instructionsstored therein causing one or more computers to perform the method ofclaim 1.